#############################################
#############################################
### REPLICATION MATERIAL                  ###  
### Gender Quotas and Political Knowledge ###
### Robustness Checks (appendix)          ###
#############################################
#############################################


# Libraries needed
library(tidyverse)
library(stargazer)
library(haven)
library(lme4)
library(rstanarm)
library(tidybayes)
library(BayesPostEst)
library(ggplot2)
library(ggridges)
library(scales)
library(rsample)
library(magrittr)
library(dplyr)
library(purrr)
library(forcats)
library(tidyr)
library(modelr)
library(ggdist)
library(cowplot)
library(rstanarm)
library(RColorBrewer)
library(gganimate)
library(wrangle)
library(ggpubr)

theme_set(theme_tidybayes() + panel_border())

# Set seed and stan options (to make it run faster)
set.seed(123)
options(mc.cores = parallel::detectCores())


#read RDS data
analysis.2 <- readRDS("data/analysis/analysis_2_df.rds")

#MODIFY FULL DATASET
analysis.int <- analysis.2 %>%
  mutate(female = recode(gender,
                         "1" = 0,
                         "0" = 1)) %>%
  select(ID, country, name, year, know_q_num, know_correct, 
         female, quota, age, education, GDP_growth, counter) %>%
  mutate(age = cut(age, 5,labels = FALSE)) %>%
  mutate(education = cut(education, 5, labels = FALSE)) %>%
  mutate(treated = ifelse(name == "Belgium" | name == "France" | name == "Greece" |
                            name == "Ireland" | name == "Portugal" |
                            name == "Spain" , 1, 0)) %>%
  write_rds("data/analysis/analysis_int.rds") %>%
  glimpse()


# Binomial version of the dataset - full sample
analysis.int.binomial <- analysis.int %>%
  group_by(female, quota, age, treated, counter, education, 
           GDP_growth, country, year, know_q_num) %>%
  summarize(n_correct = sum(know_correct),
            n_incorrect = n() - sum(know_correct)) %>%
  write_rds("data/analysis/analysis_int_binomial.rds") %>%
  glimpse()


############################################################
############################################################
### ROBUSTNESS CHECK: DIFFERENCE-IN-DIFFERENCES ANALYSIS ###
############################################################
############################################################

####################################
### MODEL AND RESULTS TABLE (A5) ###
####################################

# Model
fit.did <- stan_glmer(cbind(n_correct, n_incorrect) ~ 
                        female*quota*age +
                        female*age*treated +
                        education + 
                        GDP_growth + 
                        (1 | country) +
                        (1 | year) +
                        (1 | know_q_num) + 
                        (1 + female | country:know_q_num), 
                      data = analysis.int.binomial,
                      family=binomial(link="logit"),
                      iter = 6000,
                      control = list(adapt_delta = 0.99))

#######################
### Summary - TABLE A5
#######################
summary(fit.did)

print(fit.did, digits = 5)
mcmcReg(fit.did, format = "latex")
mcmcTab(fit.did)

#######################
### Convergence checks
#######################

# Rhat - convergence ok
dev.copy(png,'rhat_fit_did.png')
plot(fit.did, "rhat")
dev.off()

# Prior summary
prior_summary(fit.did)

#################################
### PREDICTED PROBS DATASETS  ###
#################################

# Dataframe with predicted probabilities (no random effects)

pred_fit.did_noRE <- crossing(country = unique(analysis.int.binomial$country), 
                              quota = c(0,1), 
                              female = c(0,1),
                              treated = c(0,1),
                              know_q_num = c(1:10),
                              age = c(1:5),
                              education = mean(analysis.int.binomial$education, 
                                               na.rm = T),
                              GDP_growth = mean(analysis.int.binomial$GDP_growth, 
                                                na.rm = T)) %>%
  add_epred_draws(fit.did, 
                  scale = "response", 
                  re_formula = NA) %>% 
  write_rds("data/analysis/pred_fit_did_noRE.rds")

## Dataframe with draws
draws <- fit.did %>%
  spread_draws(b[term,group]) %>%
  filter(str_detect(group, "country:know_q_num:")) %>%
  separate(group, c("country_2", "q_num_2", "country", "q_num"), ":") %>%
  write_rds("/data/analysis/draws_fit_did.rds")

# Dataframe with predicted probabilities (with random effects)

pred_fit.did_c_RE <- crossing(country = unique(analysis.int.binomial$country), 
                              quota = c(0,1), 
                              female = c(0,1),
                              year = 2005,
                              treated = c(0,1),
                              know_q_num = c(1:10),
                              age = c(1:5),
                              education = mean(analysis.int.binomial$education, 
                                               na.rm = T),
                              GDP_growth = mean(analysis.int.binomial$GDP_growth, 
                                                na.rm = T)) %>%
  add_epred_draws(fit.did, 
                  scale = "response", 
                  re_formula = ~ (1 + female | country:know_q_num)) %>%
  write_rds("/data/analysis/pred_fit_did_c_RE.rds")

#############
#############
### PLOTS ###
#############
#############

###########################################
### Figure A4: Predicted Probabilites   ###
###########################################

plot.did <- readRDS("data/analysis/pred_fit_did_noRE.rds")

plot.did$female <- factor(plot.did$female,
                          labels = c("Male", "Female"))
plot.did$quota <- factor(plot.did$quota, 
                         labels = c("No Quotas", "Quotas"))
plot.did$age <- factor(plot.did$age, 
                       labels = c("15-31", "32-48", "49-65", "66-82", "83-99"))

plot.did %>%
  ggplot(aes(x = female, y = .epred, fill = female, shape = female)) +
  stat_halfeye() +
  facet_grid(quota ~ .) +
  coord_flip() +
  scale_shape_manual(name = "Gender",
                     labels = c("Male", "Female"),
                     values = c(15, 17)) +
  scale_fill_brewer(name = "Gender",
                    labels = c("Male", "Female")) +
  labs(x = "Gender", y = "Predicted Probabilities")
ggsave("plots/appendix/fig-DID1.png")

###########################################################
### Figure A5: Predicted Probabilites - Quotas and Age  ###
###########################################################

plot.did <- readRDS("data/analysis/pred_fit_did_noRE.rds")

plot.did$female <- factor(plot.did$female,
                          labels = c("Male", "Female"))
plot.did$quota <- factor(plot.did$quota, 
                         labels = c("No Quotas", "Quotas"))
plot.did$age <- factor(plot.did$age, 
                       labels = c("15-31", "32-48", "49-65", "66-82", "83-99"))

plot.did %>%
  ggplot(aes(x = quota, y = .epred, fill = female, shape = female)) +
  stat_halfeye() +
  facet_grid(age ~ .) +
  coord_flip() +
  scale_shape_manual(name = "Gender",
                     labels = c("Male", "Female"),
                     values = c(15, 17)) +
  scale_fill_brewer(name = "Gender",
                    labels = c("Male", "Female")) +
  labs(x = "Quotas", y = "Predicted Probabilities")
ggsave("plots/appendix/fig_DID2.png")

###############################################################
### Figure A6: Predicted Probabilites - Quotas and Age (v2) ###
###############################################################

plot.did <- readRDS("data/analysis/pred_fit_did_noRE.rds")

plot.did$female <- factor(plot.did$female,
                          labels = c("Male", "Female"))
plot.did$quota <- factor(plot.did$quota, 
                         labels = c("No Quotas", "Quotas"))
plot.did$age <- factor(plot.did$age, 
                       labels = c("15-31", "32-48", "49-65", "66-82", "83-99"))

plot.did %>%
  ggplot(aes(x = age, y = .epred, fill = female, shape = female)) +
  stat_halfeye() +
  facet_grid(quota ~ .) +
  coord_flip() +
  scale_shape_manual(name = "Gender",
                     labels = c("Male", "Female"),
                     values = c(15, 17)) +
  scale_fill_brewer(name = "Gender",
                    labels = c("Male", "Female")) +
  labs(x = "Age Groups", y = "Predicted Probabilities")

ggsave("plots/appendix/fig_DID3.png")

###############################################
### Figure A7: Percent Change - Age Group   ###
###############################################

plot.did.2 <- readRDS("data/analysis/pred_fit_did_c_RE.rds")

percent_change <- plot.did.2 %>%
  ungroup() %>% 
  mutate(quota = case_when(quota == 1 ~ "quota",
                           quota == 0 ~ "noquota"),
         female = case_when(female == 0 ~ "man",
                            female == 1 ~ "woman")) %>%
  select(-.row) %>% 
  pivot_wider(names_from = c(quota, female), values_from = .epred) %>%
  mutate(second_difference = (quota_woman - quota_man) - (noquota_woman - noquota_man)) %>%
  mutate(difference_noquota = (noquota_woman - noquota_man)) %>%
  mutate(percent_change = second_difference/difference_noquota) %>%
  mutate(quotas = ifelse(country == 1 | country == 2 | country == 8 | country == 10 | country == 11 | country == 12, 1, 0)) %>%
  glimpse()

percent_change$country <- factor(percent_change$country,
                                 labels = c("France", "Belgium", 
                                            "Netherlands", "Germany", "Italy",
                                            "Luxembourg", "Denmark", "Ireland",
                                            "United Kingdom", "Greece", "Spain",
                                            "Portugal"))

percent_change$age <- factor(percent_change$age, 
                             labels = c("15-31", "32-48", "49-65", "66-82", "83-99"))

#colors
nb.cols <- 10
mycolors <- colorRampPalette(brewer.pal(8, "Dark2"))(nb.cols)

percent_change %>% 
  ggplot(aes(x = percent_change, 
             y = reorder(age, percent_change),
             color = age,
             shape = age)) + 
  stat_pointinterval() +
  coord_cartesian(xlim=c(-0.6,0.3)) +
  scale_shape_manual(name = "Age Groups",
                     labels = c("15-31", "32-48", "49-65", "66-82", "83-99"),
                     values = c(4, 15, 16, 17, 18)) +
  scale_color_manual(name = "Age Groups",
                     values = mycolors) +
  geom_vline(xintercept = 0) +
  labs(x = "Percent Change",
       y = "Age")

ggsave("plots/appendix/fig_DID5-new.png")

######################################## 
######################################## 
### PARALLEL TRENDS ASSUMPTION PLOTS ###
########################################
########################################

####################################
### Figure A8: Gender Gap Trends ###
####################################

analysis.2 <- readRDS("data/analysis/analysis_2_df.rds")

gender.know <- analysis.2 %>%
  mutate(female = recode(gender,
                         "1" = 0,
                         "0" = 1)) %>%
  spread(know_q_num,know_correct) %>%
  rename(q1 = "1", q2 = "2", q3 = "3", q4 = "4", q5 = "5",
         q6 = "6", q7 = "7", q8 = "8", q9 = "9", q10 = "10") %>%
  mutate(correct = 0) %>%
  mutate(num_question = 0) %>%
  mutate(num_question = ifelse(year == 1992, 10, num_question)) %>%
  mutate(correct = ifelse(year == 1992, q1 + q2 + q3 + q4 + q5 + q6 +
                            q7 + q8 + q9 + q10 , correct)) %>%
  mutate(num_question = ifelse(year == 1993, 4, num_question)) %>%
  mutate(correct = ifelse(year == 1993, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 1994, 4, num_question)) %>%
  mutate(correct = ifelse(year == 1994, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 1996, 10, num_question)) %>%
  mutate(correct = ifelse(year == 1996, q1 + q2 + q3 + q4 + q5 + q6 +
                            q7 + q8 + q9 + q10 , correct)) %>%
  mutate(num_question = ifelse(year == 1998, 7, num_question)) %>%
  mutate(correct = ifelse(year == 1998, q1 + q2 + q3 + q4 + q5 + q6 +
                            q7, correct)) %>%
  mutate(num_question = ifelse(year == 1999, 4, num_question)) %>%
  mutate(correct = ifelse(year == 1999, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 2002, 5, num_question)) %>%
  mutate(correct = ifelse(year == 2002, q1 + q2 + q3 + q4 + 
                            q5, correct)) %>%
  mutate(num_question = ifelse(year == 2003, 10, num_question)) %>%
  mutate(correct = ifelse(year == 2003, q1 + q2 + q3 + q4 + q5 + q6 +
                            q7 + q8 + q9 + q10, correct)) %>%
  mutate(num_question = ifelse(year == 2004, 6, num_question)) %>%
  mutate(correct = ifelse(year == 2004, q1 + q2 + q3 + q4 +
                            q5, correct)) %>%
  mutate(num_question = ifelse(year == 2005, 4, num_question)) %>%
  mutate(correct = ifelse(year == 2005, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 2006, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2006, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2007, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2007, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2008, 4, num_question)) %>%
  mutate(correct = ifelse(year == 2008, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 2009, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2009, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2010, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2010, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2011, 2, num_question)) %>%
  mutate(correct = ifelse(year == 2011, q1 + q2, correct)) %>%
  mutate(num_question = ifelse(year == 2012, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2012, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2013, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2013, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2014, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2014, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2015, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2015, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2016, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2016, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2017, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2017, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2018, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2018, q1 + q2 + q3, correct)) %>%
  mutate(fract_correct = correct/num_question) %>%
  select(ID, name, year, female, num_question, fract_correct, correct) %>%
  write_rds("data/analysis/gender_know.rds") %>%
  glimpse()

point_estimates <- gender.know %>%
  group_by(year) %>%
  ungroup() %>%
  group_by(name, year, female) %>%
  summarize(ave = mean(fract_correct)) %>% 
  pivot_wider(names_from = female, values_from = ave) %>%
  mutate(gender_gap = `0` - `1`) %>%
  select(-`0`, -`1`) 

bootstrap_function <- function() {
  gender.know %>%
    group_by(year) %>%
    sample_frac(1, replace = TRUE) %>%
    ungroup() %>%
    group_by(name, year, female) %>%
    summarize(ave = mean(fract_correct)) %>% 
    pivot_wider(names_from = female, values_from = ave) %>%
    mutate(gender_gap = `0` - `1`) %>%
    select(-`0`, -`1`) 
}

bs_iter <- 40
bs_estimates <- 1:bs_iter %>%
  map(~ bootstrap_function()) %>%
  imap(~ mutate(.x, bs_iter = .y)) %>%
  bind_rows() %>%
  glimpse()

point_estimates %>%
  ggplot(aes(x = year, y = gender_gap)) +
  geom_line(data = bs_estimates, 
            aes(x = year, y = gender_gap, group = bs_iter), 
            color = "black", alpha = 0.1) + 
  geom_line() + 
  geom_point() + 
  facet_wrap(. ~ name) +
  theme(axis.text.x=element_text(angle=60, hjust=1)) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(x = "Year",
       y = "Gender Gap",
       caption = "Gender Gap = Fraction Correct for Men - Fraction Correct for Women") +
  geom_vline(data=filter(point_estimates, name == "Belgium"), 
             aes(xintercept = 1999), colour="red") +
  geom_vline(data=filter(point_estimates, name == "France"), 
             aes(xintercept = 2000), colour="red") +
  geom_vline(data=filter(point_estimates, name == "Greece"), 
             aes(xintercept = 2012), colour="red") +
  geom_vline(data=filter(point_estimates, name == "Portugal"), 
             aes(xintercept = 2006), colour="red") +
  geom_vline(data=filter(point_estimates, name == "Spain"), 
             aes(xintercept = 2007), colour="red") +
  geom_vline(data=filter(point_estimates, name == "Ireland"), 
             aes(xintercept = 2012), colour="red")



##################################################
### Figure A9: Parallel Trend Plots by Country ###
##################################################

## Gender Gap for countries
analysis.2 <- readRDS("data/analysis/analysis_2_df.rds")

gender.know <- analysis.2 %>%
  mutate(female = recode(gender,
                         "1" = 0,
                         "0" = 1)) %>%
  spread(know_q_num,know_correct) %>%
  rename(q1 = "1", q2 = "2", q3 = "3", q4 = "4", q5 = "5",
         q6 = "6", q7 = "7", q8 = "8", q9 = "9", q10 = "10") %>%
  mutate(correct = 0) %>%
  mutate(num_question = 0) %>%
  mutate(num_question = ifelse(year == 1992, 10, num_question)) %>%
  mutate(correct = ifelse(year == 1992, q1 + q2 + q3 + q4 + q5 + q6 +
                            q7 + q8 + q9 + q10 , correct)) %>%
  mutate(num_question = ifelse(year == 1993, 4, num_question)) %>%
  mutate(correct = ifelse(year == 1993, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 1994, 4, num_question)) %>%
  mutate(correct = ifelse(year == 1994, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 1996, 10, num_question)) %>%
  mutate(correct = ifelse(year == 1996, q1 + q2 + q3 + q4 + q5 + q6 +
                            q7 + q8 + q9 + q10 , correct)) %>%
  mutate(num_question = ifelse(year == 1998, 7, num_question)) %>%
  mutate(correct = ifelse(year == 1998, q1 + q2 + q3 + q4 + q5 + q6 +
                            q7, correct)) %>%
  mutate(num_question = ifelse(year == 1999, 4, num_question)) %>%
  mutate(correct = ifelse(year == 1999, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 2002, 5, num_question)) %>%
  mutate(correct = ifelse(year == 2002, q1 + q2 + q3 + q4 + 
                            q5, correct)) %>%
  mutate(num_question = ifelse(year == 2003, 10, num_question)) %>%
  mutate(correct = ifelse(year == 2003, q1 + q2 + q3 + q4 + q5 + q6 +
                            q7 + q8 + q9 + q10, correct)) %>%
  mutate(num_question = ifelse(year == 2004, 6, num_question)) %>%
  mutate(correct = ifelse(year == 2004, q1 + q2 + q3 + q4 +
                            q5, correct)) %>%
  mutate(num_question = ifelse(year == 2005, 4, num_question)) %>%
  mutate(correct = ifelse(year == 2005, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 2006, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2006, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2007, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2007, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2008, 4, num_question)) %>%
  mutate(correct = ifelse(year == 2008, q1 + q2 + q3 + q4, correct)) %>%
  mutate(num_question = ifelse(year == 2009, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2009, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2010, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2010, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2011, 2, num_question)) %>%
  mutate(correct = ifelse(year == 2011, q1 + q2, correct)) %>%
  mutate(num_question = ifelse(year == 2012, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2012, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2013, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2013, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2014, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2014, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2015, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2015, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2016, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2016, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2017, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2017, q1 + q2 + q3, correct)) %>%
  mutate(num_question = ifelse(year == 2018, 3, num_question)) %>%
  mutate(correct = ifelse(year == 2018, q1 + q2 + q3, correct)) %>%
  mutate(fract_correct = correct/num_question) %>%
  select(ID, name, year, female, num_question, fract_correct, correct) %>%
  write_rds("data/analysis/gender_know.rds") %>%
  glimpse()

## Point estimates for countries
point_estimates <- gender.know %>%
  group_by(year) %>%
  ungroup() %>%
  group_by(name, year, female) %>%
  summarize(ave = mean(fract_correct)) %>% 
  pivot_wider(names_from = female, values_from = ave) %>%
  mutate(gender_gap = `0` - `1`) %>%
  select(-`0`, -`1`) %>%
  mutate(never_treat = ifelse(name == "Denmark" | name == "Germany" | name == "Italy" |
                                name == "Luxembourg" | name == "Netherlands" |
                                name == "United Kingdom" , 1, 0)) %>%
  mutate(treat = ifelse(never_treat == 0, 1, 0)) %>%
  glimpse()


## Belgium (1999) ##

belgium <- point_estimates %>%
  filter(name == "Belgium") %>%
  glimpse()

never_treat <- point_estimates %>%
  filter(never_treat == 1) %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

late_adopt_1 <- point_estimates %>%
  filter(treat == 1 & name != "Belgium") %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

p_belgium <- ggplot(data = point_estimates, aes(x = year)) +
  geom_line(data = belgium, aes(x = year, y = gender_gap),
            color = "#7fcdbb",
            size = 1) +
  geom_line(data = never_treat, aes(x = year, y = gender_gap_ave),
            linetype = "dashed", 
            color = "#2c7fb8",
            size = 1) +
  geom_line(data = late_adopt_1, aes(x = year, y = gender_gap_ave),
            linetype = "dotted",
            color = "#253494",
            size = 1.2) +
  geom_vline(xintercept = 1999) +
  labs(x = "Year",
       y = "Gender Gap") +
  ggtitle("Belgium") +
  theme(plot.title = element_text(hjust = 0.5))
p_belgium

## France (2000)##

france <- point_estimates %>%
  filter(name == "France") %>%
  glimpse()

never_treat <- point_estimates %>%
  filter(never_treat == 1) %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

late_adopt_2 <- point_estimates %>%
  filter(treat == 1 & name != "Belgium" & name != "France") %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

early_adopt_1 <- point_estimates %>%
  filter(treat == 1 & name == "Belgium") %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

p_france <- ggplot(data = point_estimates, aes(x = year)) +
  geom_line(data = france, aes(x = year, y = gender_gap),
            color = "#7fcdbb",
            size = 1) +
  geom_line(data = never_treat, aes(x = year, y = gender_gap_ave),
            linetype = "dashed", 
            color = "#2c7fb8",
            size = 1) +
  geom_line(data = late_adopt_2, aes(x = year, y = gender_gap_ave),
            linetype = "dotted",
            color = "#253494",
            size = 1.2) +
  geom_line(data = early_adopt_1, aes(x = year, y = gender_gap_ave),
            linetype = "dotdash",
            color = "#41b6c4",
            size = 1) +
  geom_vline(xintercept = 2000) +
  labs(x = "Year",
       y = "Gender Gap") +
  ggtitle("France") +
  theme(plot.title = element_text(hjust = 0.5))
p_france

## Portugal (2006) ##

portugal <- point_estimates %>%
  filter(name == "Portugal") %>%
  glimpse()

never_treat <- point_estimates %>%
  filter(never_treat == 1) %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

late_adopt_3 <- point_estimates %>%
  filter(treat == 1 & name != "Belgium" & name != "France" & name != "Portugal") %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

early_adopt_2 <- point_estimates %>%
  filter(treat == 1 & name == "Belgium" | name == "France") %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

p_portugal <- ggplot(data = point_estimates, aes(x = year)) +
  geom_line(data = portugal, aes(x = year, y = gender_gap),
            color = "#7fcdbb",
            size = 1) +
  geom_line(data = never_treat, aes(x = year, y = gender_gap_ave),
            linetype = "dashed",
            color = "#2c7fb8",
            size = 1) +
  geom_line(data = late_adopt_3, aes(x = year, y = gender_gap_ave),
            linetype = "dotted",
            color = "#253494",
            size = 1.2) +
  geom_line(data = early_adopt_2, aes(x = year, y = gender_gap_ave),
            linetype = "dotdash",
            color = "#41b6c4",
            size = 1) +
  geom_vline(xintercept = 2006) +
  labs(x = "Year",
       y = "Gender Gap") +
  ggtitle("Portugal") +
  theme(plot.title = element_text(hjust = 0.5))
p_portugal

## Spain (2007) ##

spain <- point_estimates %>%
  filter(name == "Spain") %>%
  glimpse()

never_treat <- point_estimates %>%
  filter(never_treat == 1) %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

late_adopt_4 <- point_estimates %>%
  filter(treat == 1 & name != "Belgium" & name != "France" & name != "Portugal" 
         & name != "Spain" ) %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

early_adopt_3 <- point_estimates %>%
  filter(treat == 1 & name == "Belgium" | name == "France" | name == "Portugal") %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

p_spain <- ggplot(data = point_estimates, aes(x = year)) +
  geom_line(data = spain, aes(x = year, y = gender_gap),
            color = "#7fcdbb",
            size = 1) +
  geom_line(data = never_treat, aes(x = year, y = gender_gap_ave),
            linetype = "dashed", 
            color = "#2c7fb8",
            size = 1) +
  geom_line(data = late_adopt_4, aes(x = year, y = gender_gap_ave),
            linetype = "dotted",
            color = "#253494",
            size = 1.2) +
  geom_line(data = early_adopt_3, aes(x = year, y = gender_gap_ave),
            linetype = "dotdash",
            color = "#41b6c4",
            size = 1) +
  geom_vline(xintercept = 2007) +
  labs(x = "Year",
       y = "Gender Gap") +
  ggtitle("Spain") +
  theme(plot.title = element_text(hjust = 0.5))
p_spain

## Greece (2012) ##

greece <- point_estimates %>%
  filter(name == "Greece") %>%
  glimpse()

never_treat <- point_estimates %>%
  filter(never_treat == 1) %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

late_adopt_5 <- point_estimates %>%
  filter(treat == 1 & name != "Belgium" & name != "France" & name != "Portugal" 
         & name != "Spain" & name != "Greece") %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

early_adopt_4 <- point_estimates %>%
  filter(treat == 1 & name == "Belgium" | name == "France" | name == "Portugal"
         | name == "Spain") %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

p_greece <- ggplot(data = point_estimates, aes(x = year)) +
  geom_line(data = greece, aes(x = year, y = gender_gap),
            color = "#7fcdbb",
            size = 1) +
  geom_line(data = never_treat, aes(x = year, y = gender_gap_ave),
            linetype = "dashed", 
            color = "#2c7fb8",
            size = 1) +
  geom_line(data = late_adopt_5, aes(x = year, y = gender_gap_ave),
            linetype = "dotted",
            color = "#253494",
            size = 1.2) +
  geom_line(data = early_adopt_4, aes(x = year, y = gender_gap_ave),
            linetype = "dotdash",
            color = "#41b6c4",
            size = 1) +
  geom_vline(xintercept = 2012) +
  labs(x = "Year",
       y = "Gender Gap") +
  ggtitle("Greece") +
  theme(plot.title = element_text(hjust = 0.5))
p_greece

## Ireland (2012) ##

ireland <- point_estimates %>%
  filter(name == "Ireland") %>%
  glimpse()

never_treat <- point_estimates %>%
  filter(never_treat == 1) %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

early_adopt_5 <- point_estimates %>%
  filter(treat == 1 & name == "Belgium" | name == "France" | name == "Portugal"
         | name == "Spain" | name == "Greece") %>%
  group_by(year) %>%
  summarize(gender_gap_ave = mean(gender_gap)) %>%
  glimpse()

p_ireland <- ggplot(data = point_estimates, aes(x = year)) +
  geom_line(data = ireland, aes(x = year, y = gender_gap),
            color = "#7fcdbb",
            size = 1) +
  geom_line(data = never_treat, aes(x = year, y = gender_gap_ave),
            linetype = "dashed", 
            color = "#2c7fb8",
            size = 1) +
  geom_line(data = early_adopt_4, aes(x = year, y = gender_gap_ave),
            linetype = "dotdash",
            color = "#41b6c4",
            size = 1) +
  geom_vline(xintercept = 2012) +
  labs(x = "Year",
       y = "Gender Gap") +
  ggtitle("Ireland") +
  theme(plot.title = element_text(hjust = 0.5))
p_ireland

## Commbine plots into 1 ##
par_trends <- ggarrange(p_belgium, p_france, p_portugal, 
                        p_spain, p_greece, p_ireland,
                        ncol = 2, nrow = 3)
par_trends

ggsave("plots/main/trends.png", 
       width = 6.5, 
       height = 7, 
       units = "in")

annotate_figure(par_trends,
                top = text_grob("Assessing Parallel Trends", 
                                color = "black", 
                                face = "bold", 
                                size = 14),
                bottom = text_grob("Legend \n Solid: Treated \n Dashed: 
                                   Never Treated \n Dotted: Late Adopters 
                                   \n Dotdashed: Early Adopters ", 
                                   color = "black",
                                   hjust = 1, x = 1, 
                                   face = "italic", 
                                   size = 5))









